%function [speed Svec]=compute_spreading(lon,lat,plates); 
clear;
load ../data/generate_profiles/collect_all.mat
load morvel.mat P wc ws;
ro=6.378*10^6*10^2; %radius of the earth in cm

%for each bathymetry profile, identify the plate pairing and compute
%velocities
for i=1:length(All.S);
    lon=All.loc(i,1);
    lat=All.loc(i,2);
    k=All.region(i);
    P=plate_pair(k,ws); 

    for ct=1:2,
        [P(ct).v P(ct).s P(ct).theta]=euler_convert(lat,lon,P(ct).ws);
    end;

    All.velocity(i,:)=[P(1).v-P(2).v]; %mm/year
    All.speed(i)=sqrt(sum(All.velocity(i,:).^2)); %mm/year
end;

figure(1); clf; hold on;
set(gca,'FontSize',16);
box on;
plot(All.S,'LineWidth',1.5);
plot(All.speed/20,'--','LineWidth',1.5);
xlabel('Segment number');
ylabel('Spreading half-rate (cm/yr)');
legend('B-M method','MORVEL solution');

figure(2); clf; hold on;
box on;
set(gca,'FontSize',16);
colormap hsv; 
S1=All.S;
S2=All.speed/20;
plot([1.5 8],[1.5 8],'k--','linewidth',1.5); 

pl=find(All.side==1);
h=scatter(S2(pl),S1(pl),50,All.region(pl),'o'); 
h.MarkerFaceColor='Flat';

pl=find(All.side==2);
h=scatter(S2(pl),S1(pl),50,All.region(pl),'o'); 
h.MarkerFaceColor='Flat';

h=xlabel('half spreading rate (cm/year MORVEL)'); 
h=ylabel('half spreading rate (cm/year this study)');
axis tight;

colormap jet;
temp=colormap;
for ct=1:3,
    c(:,ct)=interp1(linspace(0,17,length(temp)),temp(:,ct),[1:17]-0.5);
end;
set(gca,'clim',[0 17]);
colormap(c);
h=colorbar;
h.Ticks=0.5:16;
h.TickLabels={1:16};
h.Label.String='region number';

